Spectral signatures of excess-proton waiting and transfer-path dynamics in aqueous hydrochloric acid solutions

The theoretical basis for linking spectral signatures of hydrated excess protons with microscopic proton-transfer mechanisms has so far relied on normal-mode analysis. We introduce trajectory-decomposition techniques to analyze the excess-proton dynamics in ab initio molecular-dynamics simulations of aqueous hydrochloric-acid solutions beyond the normal-mode scenario. We show that the actual proton transfer between two water molecules involves for relatively large water-water separations crossing of a free-energy barrier and thus is not a normal mode, rather it is characterized by two non-vibrational time scales: Firstly, the broadly distributed waiting time for transfer to occur with a mean value of 200–300 fs, which leads to a broad and weak shoulder in the absorption spectrum around 100 cm−1, consistent with our experimental THz spectra. Secondly, the mean duration of a transfer event of about 14 fs, which produces a rather well-defined spectral contribution around 1200 cm−1 and agrees in location and width with previous experimental mid-infrared spectra.

T he motion of excess protons in aqueous solution is fundamental for many biological and chemical processes. The excess-proton diffusivity is significantly higher compared to other monovalent cations in water 1,2 , since the excess proton exchanges its identity with water hydrogens during the diffusion process 3,4 . Grotthus hypothesized a similar process over two centuries ago 5,6 , but a detailed understanding of the protontransfer dynamics between water or other molecules remains difficult to date due to the multitude of time scales involved and the only indirect experimental evidence.
Infrared (IR) spectroscopy in the THz and mid-IR regimes is a powerful tool to explore the ultrafast dynamics of water and aqueous ion solutions. For example, the prominent absorption peak around 200 cm −1 of bulk water is dominated by firstsolvation-shell dynamics, whereas motion involving the second solvation shell contributes most significantly below 80 cm −1 (2.4 THz) 7,8 . Furthermore, so-called 'rattling' modes for strongly hydrated ions lead to characteristic absorption features, while for weakly hydrated ions vibrationally induced charge fluctuations are dominant 9,10 , as suggested by dissecting simulation spectra into contributions from different solvation shells 11,12 .
IR spectroscopy has proven particularly useful for the study of the ultrafast dynamics of the excess proton in aqueous solution 13,14 . Due to their low pH value, aqueous hydrochloric acid (HCl) solutions are perfect model systems to study excessproton dynamics, since HCl dissociates readily in water and gives rise to a large number of highly mobile protons that are only weakly coordinated with neutralizing chloride ions and therefore behave as if added in excess. The characteristic continuum band in the IR absorption spectrum, located between the waterbending mode around 1650 cm −1 and the water-stretching mode around 3300 cm −1 , has long been known and led to the hypothesis of the Zundel state, i.e., two water molecules symmetrically sharing the excess proton 15 . This model has been challenged by a contrasting picture, the Eigen state, which is a hydronium ion caged symmetrically by three water molecules 16 . Ever since these idealized structures have been proposed, their relative stability has been controversially debated 13,[17][18][19][20][21][22][23][24][25][26][27] . It is now known that neither the idealized Zundel nor the Eigen states are realistic structural representations and that the excess proton mostly resides slightly asymmetrically shared between two water molecules, in the 'special pair' state, which geometrically can be interpreted as a 'distorted Zundel' state or a 'distorted Eigen' state 20,23,24,[26][27][28] .
While simulations can reproduce most experimental spectroscopic signatures, the understanding of the proton transfer mechanism requires model building based on and guided by simulations. It is generally accepted that proton transfer involves consecutive transitions between states that can be viewed as more Eigen-like and more Zundel-like and have fast interconversion times 24,26,29 . That the excess proton diffusion involves the crossing of free-energetic barriers follows from the experimentally known Arrhenius behavior of the excess-proton conductivity 28,[30][31][32] . Similar to the above-mentioned discussion on the relative stability of the Eigen and Zundel states, it remains debated whether the Zundel state is the transition state between two Eigen states or the opposite is the case, i.e., whether the Eigen state is the transition state between two Zundel states 23,24,27 . Theoretical models for the spectroscopic signatures of the hydrated proton motion so far relied on normal-mode calculations and have explained many aspects of experimental linear absorption [33][34][35] as well as 2D IR spectra 25,26 . However, normal modes by construction cannot deal with the thermally activated transfer of an excess proton over a free-energy barrier, since this corresponds to an unstable mode with a negative free-energy curvature along the transfer reaction coordinate 36 . It is clear that such proton-barrier-transfer events will make a sizable spectroscopic contribution, since they involve fast motion of a highly charged object over relatively large distances. From this follows that an excess-proton transition state, which corresponds to a free energy maximum and thus occurs with a small probability, nevertheless can make a dominant contribution to the spectrum, which would lead to characteristic differences between experimental spectra and normal-mode theory predictions. Indeed, it has been noted that the normal-mode spectra computed from instantaneous configurations do not explain all experimental spectral signatures associated with the excess proton in water 19,[25][26][27] , in particular of the proton-transfer dynamics 37 . It was recently shown that proton transfer in the H 5 O 2 + cationic complex gives rise to two distinct spectroscopically relevant time scales that cannot be captured by normal-mode analysis 38 .
In general, the transfer of a particle with mass m over an energy barrier with negative curvature k < 0 corresponds to an unstable mode. The dynamics of such a barrier-crossing is not characterized by a vibrational time scale, which according to a harmonic oscillator model could erroneously be written as , where U 0 is the barrier height and L is the barrier width, but rather by two other time scales, namely the mean transfer-waiting time τ TW and the mean transfer-path time τ TP . τ TW is the average waiting time before a transfer event occurs and τ TP is the average time of the actual transfer path over the energy barrier. The former one scales exponentially with the barrier height U 0 , τ TW $ expðU 0 =k B TÞ 36,39 , whereas the latter one scales logarithmically with the barrier height U 0 , τ TP $ logðU 0 =k B TÞ [40][41][42][43] . In essence, it is not clear with current theoretical methodology what the spectroscopic signature of excess-proton transfer events in aqueous solutions is and whether the continuum band stems just from vibrations in metastable states or whether transfer reactions over barriers are involved. Thus, a theoretical approach that is complementary to normal modes and can handle proton-transfer events that involve freeenergy barriers is needed.
In this study, we investigate the excess-proton dynamics in aqueous HCl solutions at ambient conditions using ab initio molecular dynamics (MD) simulations at the Born-Oppenheimer level and experimental THz/Fourier-transform infrared (THz/ FTIR) measurements. Our simulated IR difference absorption spectra compare well to our experimental data in the THz regime as well as to literature data in the mid-IR regime. By projecting the excess-proton dynamics onto the two-dimensional coordinate system spanned by the proton position along the axis connecting the two closest water oxygens d and the oxygen distance R OO 44,45 , the excess-proton trajectories and their spectral signatures are subdivided into three contributions with distinct time scales, as illustrated in Fig. 1a. The fastest time scale, τ NM , reflects vibrations when the excess proton transiently forms a solvated H 3 O + molecule which is asymmetrically solvated in a special pair. It is well captured by a normal-mode description and has been amply discussed in literature 19,25,26,37,46 . The other two spectral signatures, stemming from proton transfer events, are the focus of this study. The associated time scales τ TW and τ TP are directly obtained from our simulated excess proton trajectories using a multidimensional path analysis. While the transfer-waiting time τ TW of aqueous proton-transfer events has been studied recently 47 , the identification of both τ TW and τ TP in simulated and experimental spectra is a main result of this work. We find a mean transfer-waiting time of τ TW = 200-300 fs depending on HCl concentration, which in our experimental THz spectra shows up as a broad weak shoulder around 100 cm −1 , that is partially overlaid by the absorption due to rattling chloride anions at about 150 cm −1 9,11 . The mean transfer-path time, from simulations obtained as τ TP = 14 fs, produces a spectroscopic signature around 1200 cm −1 , which is well captured in experimental mid-IR spectra 14,20,23,24,26 . Note that in our simulations the proton transfer becomes barrier-less for small water separation and thus includes the highly anharmonic normal-mode vibration of Zundel-like configurations, which have already been analyzed theoretically 19,25,37,46 . In the THz regime our experimental difference spectra show an additional prominent peak around 300-400 cm −1 , in good agreement with our simulation data, which is demonstrated to be caused by the coupling of the excessproton motion to the relative oscillations of the two flanking water molecules in transient H 5 O 2 + complexes. Proton transfer events between water molecules are frequently followed by a few immediate back-and-forth transfer events, which is a consequence of non-Markovian effects 48 that have to do with the slowly changing solvation structure around the excess proton. Although these transfer events are therefore not always productive in the sense that they lead to large-scale diffusion of the excess proton, they nevertheless give rise to pronounced experimental spectroscopic signatures and therefore need to be included in the analysis.

RESULTS
Infrared and THz spectra of HCl solutions. Within linear spectroscopy, the energy absorption rate of incident light with frequency f = ω/(2π) is proportional to the imaginary part of the dielectric susceptibility and given by ωχ″(ω). IR power spectra are obtained from ab initio MD simulations of water (blue solid line) and HCl solutions at three concentrations between 2 and 6 M (purple to red solid lines) and are shown in Fig. 1b. The spectra are divided by the water molecular number concentration ωχ 00 c W ¼ ωχ 00 =c W . Simulation details are provided in the "Methods" section. All IR spectra show the characteristic features of pure water spectra, which are the prominent OH-stretching peak around 3300 cm −1 , the HOH-bending mode around 1650 cm −1 and librational modes in the far IR regimes between 200 and 800 cm −1 . The IR spectra of HCl solutions additionally show a broad continuum between the bending and the stretching peaks, from 2000 to 3000 cm −1 , and a broad peak at around 1200 cm −1 , both of which are commonly interpreted as to reflect the excessproton dynamics 14,23,49 . Furthermore, additional features are observed below 800 cm −1 , that are shown in Fig. 2 in comparison to our experimental THz spectra and will be discussed further below.
The simulated difference spectra in Fig. 1c (solid lines) clearly demonstrate three distinct regions (color shaded), that relate to distinct time scales of the excess-proton dynamics and will in this work be identified as transfer-waiting (TW, gray), transfer-path (TP, red) and normal-mode contributions (NM, green). We obtain rather good agreement with the experimental difference spectrum for 4 M HCl 14 , which was scaled to match the height of the simulated IR 1200 cm −1 peak, see Supplementary Fig. 5 for a comparison of different experimental data. Our simulated 4 M difference spectrum in Fig. 1c does not reproduce the local maximum of the experimental difference spectra around 1750 cm −1 , which is interpreted as the acid-bend band, i.e., a blue shift of the bending mode in H 3 O + compared to water, and also not the shape of the experimental acid- Fig. 1 Time scales of excess-proton dynamics and simulated absorption spectra. a Schematic trajectory of an excess proton that transfers between two water molecules, together with a schematic free energy profile F(d) that exhibits a barrier and is representative of a relatively large oxygen-oxygen separation R OO . Three time scales characterize the proton trajectory, the normal-mode vibrational period of the solvated transient H 3 O + , τ NM , the transferpath time, τ TP and the transfer-waiting time, τ TW , where τ TW > τ TP > τ NM . An animation is shown online https://fu-berlin.eu.vbrickrev.com/sharevideo/ df2d94a4-6e7f-499a-a256-17d73b6124e4. b Infrared (IR) absorption spectra obtained from ab initio molecular dynamics (MD) simulations of pure water (blue solid line) and hydrochloric acid (HCl) solutions at various concentrations (dark purple: 2 M, purple: 4 and red: 6 M). The spectra are divided by the water molecular number concentration c W . c Difference spectra between the three HCl spectra and the water spectrum, obtained from the spectra in b. The purple dotted line shows an experimental difference spectrum of HCl at 4 M 14 , rescaled in height to match the simulation results. d The simulated difference spectra (as shown in c) divided by the HCl concentrations c HCl . Three distinct spectral regions are shaded in different colors, that are identified with different excess-proton dynamic processes: transfer waiting (TW, gray), transfer paths (TP, red), and normal modes (NM, green). The transfer-waiting time is close to the chloride-ion (Cl − ) rattling time and the oxygen vibrational time in local H 5 O 2 + complexes that is described by the R OO coordinate. stretch signature around 3000 cm −1 14,19,20,23,50 . The reason for this disagreement is unclear, we note that the normalization of spectra when calculating difference spectra is a subtle issue, see Supplementary Note 3 for a discussion. Figure 1d shows that the three simulated HCl difference spectra divided by the HCl concentrations c HCl are nearly indistinguishable. This clearly indicates that the spectroscopic features are due to single-proton dynamics and that collective proton effects as well as proton-chloride coupling effects, which would scale non-linearly in the HCl concentration, are minor. This is an important finding and justifies our theoretical analysis of single excess-proton motion in this work.
In order to investigate the intermolecular vibrational dynamics of water, solvated protons, and chloride ions, we experimentally measure THz absorption spectra for HCl concentrations of 2 M, 4 M, 6 M, 8 M, and 10 M. For quantitative comparison with our simulation data, the experimentally measured extinction spectra are converted into energy absorption spectra using the Kramers-Kronig relation, details are described in the "Methods" section and in Supplementary Methods 1. The experimental THz/ FTIR spectra are shown in Fig. 2a in the range 0-650 cm −1 (colored broken lines) together with a literature spectrum of pure water (blue broken line) and are compared to the available simulated spectra (solid lines). Again, all experimental and simulated spectra are divided by the respective water concentration c w . One notes the good agreement between the experimental and simulation spectra below 400 cm −1 , which is noteworthy since the spectral amplitudes are not rescaled or adjusted. However, the reason of the disagreement for larger wave numbers is not clear. All spectra show a prominent peak at 200 cm −1 . Difference spectra of the experimental data with respect to the pure water spectrum are shown in Fig. 2b (broken lines) and again compared to the available simulated difference spectra (solid lines). Two peaks dominate the difference spectra, one around 150-200 cm −1 and one around 300-400 cm −1 . The experimental difference spectra scale linearly with HCl concentration, which is demonstrated in Fig. 2c, where the difference spectra are divided by the HCl concentrations c HCl . For comparison, the simulated difference spectra divided by c HCl , already presented in Fig. 1d, are averaged over the three HCl concentrations and shown as a black solid line. The linear scaling of the experimental spectra with HCl concentration reconfirms that the difference spectra are related to single-ion behavior and that collective ion effects are negligible, in agreement with previous observations 13, 14 . In essence, two different processes at 150-200 cm −1 and 300-400 cm −1 are clearly indicated by our experimental and simulated spectra and will be interpreted by our spectral trajectory-decomposition techniques. In the remainder we analyze exclusively the 6 M solution, which provides the best proton statistics.
Excess-proton trajectories and spectra. Excess protons constantly change their identity as they move through the HCl solution. Each identity change introduces a spurious discontinuity in the excess-proton trajectory, which does not actually correspond to charge transport and therefore is spectroscopically irrelevant. In order to extract continuous excess-proton trajectories from our simulations, we use a dynamic criterion as illustrated in Fig. 3a (that our extracted excess-proton trajectories are spectroscopically meaningful we will a posteriori demonstrate by comparison of spectra calculated from excess proton trajectories with spectra calculated from the complete simulation system). Each proton is assigned to its closest oxygen atom at each time step. Whenever three protons are assigned to the same oxygen, thereby forming a hydronium ion, all of them are registered as excess-proton candidates. That means, for the generation of continuous excess-proton trajectories, we do not select the hydrogen with the largest separation from the oxygen, which would lead to fast switching of the excess proton identity, the socalled 'special pair dance' of hydronium with its surrounding water molecules 27,51 . Rather, if during the simulation an excessproton candidate becomes assigned to a different oxygen and thus transfers to a neighboring water, it is selected as an excess proton for the entire time during which it was part of any hydronium ion 27 . Note that the spectral effects of the rattling of the excessproton candidates within one hydronium ion, i.e., the 'special pair dance', are in some of our calculations below included by taking into account the flanking water molecules in the calculation of spectra, but do not show significant spectral signatures. Excess protons that are coordinated with a chloride anion as the second nearest neighbor are neglected from our analysis. This does not influence our excess-proton spectra, since even for the highest acid concentration of 6 M, only 5% of all configurations are of this type, as demonstrated in Supplementary Table 3. Note, however, that the fraction of protons coordinated with chloride ions increases significantly at higher concentrations 52 . Our procedure for calculating continuous excess-proton trajectories is discussed in further detail in Supplementary Methods 2.
The excess-proton trajectories are described by the twodimensional coordinate system defined within local transient H 5 O 2 + complexes consisting of the excess proton and its two nearest water molecules, as illustrated in the right part of Fig. 3a. The coordinates are the instantaneous distance between the two Fig. 2 Experimental absorption spectra. a Experimental THz/Fourier-transform infrared (THz/FTIR) absorption spectra of HCl solutions at various concentrations (colored broken lines for 2 M to 10 M), compared to literature data of pure water 74 (blue broken line). The experimentally measured extinction spectra have been converted into energy absorption spectra using the Kramers-Kronig relation, no amplitude adjustment is used in the comparison with the ab initio spectra (solid lines). b Experimental and ab initio molecular dynamics (MD) difference spectra derived from the results given in a and plotted in the respective colors and line styles. c The experimental difference spectra (shown in b) are divided by the HCl concentration c HCl (colored broken lines) and compared to the average of the simulated difference spectra after dividing by c HCl (black solid line).
oxygen atoms R OO and the excess proton's distance from the midplane d 45 . The state for d = 0, where the excess proton is in the middle between the oxygens, will be later used to define the transition state of the proton transfer between the two flanking water molecules. Figure 3b shows an example excess-proton trajectory from our ab intio MD simulations in terms of the R OO and d coordinates. While the motion along the two coordinates is strongly correlated, as we will show later, the d trajectory shows fast oscillatory components that are much weaker in the R OO trajectory.
The power spectra of the excess-proton trajectories, averaged over all excess protons in the solution, are shown in Fig. 3c for the d coordinate as a black broken line and for the R OO coordinate as a black dotted line, in the calculations we assume a bare charge of 1 e for the excess proton (left axis). We compare with the simulated difference spectrum of the 6 M HCl solution (red solid line), which is multiplied by the water concentration c W and divided by HCl concentration c HCl and thus is normalized per excess proton (right axis). The qualitative agreement between the two spectra (black broken and red solid lines) is very good up to an overall scaling factor of roughly four, which reflects polarization enhancement due to neighboring water molecules. The good agreement indicates that the difference spectrum of an HCl solution is proportional to the spectrum of the highly IRactive excess proton in terms of its coordinate d 47 . In other words, the HCl-solution difference spectrum reports on the excessproton motion relative to the two flanking water oxygens and can therefore be used to investigate proton-transfer dynamics. In turn, the analysis of excess-proton trajectories allows to reveal the microscopic mechanism causing the signatures of HCl-solution difference spectra, which is a central validation of the trajectorydecomposition technique used in our study. In contrast, the dynamics of R OO , i.e., the vibrations of the water molecules in the H 5 O 2 + complex, black dotted line in Fig. 3c, gives rise to a single spectral feature around 400 cm −1 , which turns out to be present also in the simulated and experimental HCl-solution spectra, as will be discussed below.
To check for the effect of the two water molecules that flank the excess proton on the difference spectrum, we also calculate the IR spectrum of transient H 5 O 2 + complexes, as done previously 22,53 and presented in more detail in Supplementary Fig. 11. To construct a difference spectrum, we subtract from the H 5 O 2 + spectrum the spectrum of hydrogen-bonded water-molecule pairs obtained from the pure-water ab initio MD simulation. The resulting difference spectrum in Fig. 3c (gray solid line, right scale) is reduced by a factor of roughly two compared to the difference spectrum of the entire HCl solution (red solid line) but otherwise agrees in shape rather nicely. Compared to the spectrum of the isolated excess proton (broken line, left scale) we observe an amplification by a factor of roughly two, but no essential spectral shape change. We conclude that the flanking water molecules and in particular the 'special pair dance' with further solvating water molecules does not modify the spectrum of the excess proton in an essential way. The amplification of the complete HCl-solution difference spectrum compared to the H 5 O 2 + difference spectrum (red and gray solid lines, respectively) we rationalize by polarization enhancement effects of water molecules that solvate the H 5 O 2 + complex.
A few spectral contributions that are not included in the excess proton power spectrum (black broken line in Fig. 3c) deserve mentioning: (i) Dynamics orthogonal to the connecting axis of the oxygens are shown to be small in Supplementary Fig. 11. (ii) The chloride motion is shown below to contribute only slightly and at low frequencies to the spectrum. iii) The translation and rotation of the internal H 5 O 2 + coordinate system relative to the lab frame is in Supplementary Figs. 12 and 13 shown to only give a small spectral contribution. We thus conclude that the IR difference spectrum between HCl solutions and pure water reports very faithfully on the excess-proton dynamics, apart from an overall amplification factor.
2D excess-proton trajectory analysis. Figure 4a shows the twodimensional (2D) free energy for 6 M HCl obtained from the negative logarithm of the distribution function of the continuous excess-proton trajectories as a function of the coordinates d and R OO , a blow up of the shaded area is given in Fig. 4e. By definition, the free energy is symmetric with respect to the midplane at d = 0, which separates two global minima at R OO = 2.51 Å and d = ±0.2 Å. These minima, highlighted as gray dots in Fig. 4a, correspond to states where the excess proton is asymmetrically shared between the two flanking water molecules. The transition between these minima, i.e., the proton transfer, is, therefore, a barrier-crossing process in the two-dimensional plane spanned by R OO and d. Figure 4b shows cuts through the free energy for constant R OO along d, each fitted to a quartic expression F(d) = F d=0 (1 + γ 2 d 2 + γ 4 d 4 ) shown as black broken lines. Details are reported in Supplementary Methods 3. For negative γ 2 , two minima at d Ã ¼ ± ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Àγ 2 =ð2γ 4 Þ p are separated by a barrier at d = 0, which correspond to the optimal proton asymmetry for a given value of R OO and are determined by the parabolic function d * (R OO ), which is plotted as a black broken line in Fig. 4a, e. The cuts in The absolute free energy at d = 0 is plotted in Fig. 4c as a red solid line and compared to the barrier height relative to the R OOdependent minima at d * (R OO ) (blue solid line). The minimal absolute barrier free energy of 0.9 k B T (red line), located at R OO = 2.42 Å, defines the transition state; for R OO = 2.51 Å, for which the most probable excess-proton state is obtained, the barrier has a moderate absolute height of 1.8 k B T, suggesting that proton transfer is not excluded for this value of R OO . Note that for R OO < 2.39 Å the relative barrier height vanishes and thus a symmetrically shared excess proton is most likely.
Next, to decompose the excess proton trajectories into segments where the excess proton moves around the local free energy minima and where a transfer across the midplane happens, transfer-path start and end points need to be defined. For this we use the most likely proton location d * (R OO ) (black broken lines in Fig. 4a, e, g). The start of a transfer path is thus defined as the last crossing of d * (R OO ) on one side of the midplane at d = 0 and the end of a transfer path as the first crossing of d * (R OO ) on the other side of the midplane at d = 0. The transfer paths are slightly extended forward and backward in time to the points where the velocity along d vanishes, the socalled turning points, in order to be consistent with the analytical theory presented in 38 . An example trajectory is shown in Fig. 4e (thin black line) in the (d, R OO ) plane, the corresponding timedependent position d is given in Fig. 4g, where the transfer path is highlighted in blue. Many transfer attempts are unsuccessful and lead to incomplete transfer paths, where the excess proton crosses the mid-plane d = 0 but does not reach to the minimal free energy state d * (R OO ) on the other side. Two incomplete transfer paths are shown in green in Fig. 4g and for consistency are also extended to their turning points. Transfer-path-time distributions are given in Fig. 4f; with our definitions used, incomplete transfer paths turn out to be slightly faster. In total, there are about as many incomplete (n = 14877) as complete transfer paths (n = 14357), meaning that about half of all excess protons reaching the midplane d = 0 actually transfer from one water molecule to the other. The main peaks in the distributions are fitted by the Erlang distribution 54-56 with the mean transfer-path time defined by τ TP , shown as black solid lines in Fig. 4f, the fit parameters are given in the legend. The distribution of transition states in Fig. 4d, i.e., the R OO position at which complete transfer paths cross the midplane at d = 0, is rather broad and peaks slightly below R OO = 2.42 Å, the most probable excess-proton position at d = 0. Most paths, in fact, 77%, cross for R OO > 2.39 Å, i.e., for values of R OO where a barrier along the d coordinate is present. This means that the dominant mechanism for proton transfer is not one where the proton waits until the oxygen-oxygen separation R OO reaches small values so that the remaining barrier along d is small or absent. Rather, protons cross the d = 0 midplane for a broad distribution of oxygen-oxygen separations R OO and by doing so overcome substantial free-energy barriers. This reverberates that a normal-mode analysis cannot account for all aspects of proton transfer in HCl solutions.
Spectral signatures of proton transfer. In order to dissect the excess-proton spectrum in Fig. 3c (black broken line) into contributions that have to do with proton-transfer events and those that do not, the excess-proton trajectories d(t) are decomposed into three parts according to To illustrate this decomposition, Fig. 5a shows part of an example excess-proton trajectory, d(t) (black line), together with the most likely excess-proton positions d * [R OO (t)] (thin gray lines); the deviations between the black and gray lines visualize excessproton motion relative to the oxygen it is bound to. We define the transfer-waiting contribution to the excess-proton trajectory as d TW (t) ≡ d * [R OO (t)] projected onto the closer branch of d * [R OO (t)], shown as a blue solid line in Fig. 5b. Thereby, d TW (t) reflects the proton transfer jumps and also contains the water motion. The transfer-path contribution d TP (t) in Fig. 5c (red solid line) is defined as d TP (t) = d(t) − d TW (t) during complete and incomplete transfer paths (as defined in Fig. 4e, g) and is zero elsewhere, it describes the excess-proton motion during transfer processes. Finally, by subtracting d TW (t) and d TP (t) from d(t), we are left with the oscillations around d * [R OO (t)] when the excess proton is not undergoing a transfer, which constitutes the normal-mode contribution d NM (t) in Fig. 5d (green solid line). Different or more detailed excess-proton trajectory decompositions are certainly conceivable, the usefulness of the present scheme follows from its spectral decomposition properties. In Fig. 5e the excess-proton spectrum, (black solid line) is decomposed as The power spectra of the transfer-waiting, χ 00 TW (blue line), and transfer-path contributions, χ 00 TP (red line), are computed from the d TW (t) and d TP (t) trajectories using the Wiener-Kintchine theorem (see "Methods" section for details). All crosscorrelation contributions are included in the normal-mode contribution, χ 00 NM (green line). The normal-mode spectrum ωχ 00 NM in Fig. 5e accounts for the continuum band located between 2000 and 3000 cm −1 , it is in fact amenable to normal-mode analysis 19,25,27 but by construction does not include the proton-transfer dynamics. The range of the dominant normal-mode time scales included in ωχ 00 NM , τ NM = 11-17 fs, follows from the spectral width of the continuum band, taken to be f = 2000-3000 cm −1 in Fig. 5e, via τ NM = 1/f. The transfer-path spectrum ωχ 00 TP in Fig. 5e shows a pronounced peak around 1200 cm −1 . An analytical model calculation shows that the peak in the transfer-path spectrum is related to the mean transfer-path time as f TP = 1/(2τ TP ) 38 . Taking the results from the fits in Fig. 4f, yielding τ c TP = 14.1 fs for complete and τ i TP = 12.7 fs for incomplete transfer paths, we predict f c TP ¼ 1170 cm À1 and f i TP ¼ 1300 cm À1 , indicated in Fig. 5e as vertical lines and which bracket the transfer-path peak very nicely.
The transfer-waiting spectrum ωχ 00 TW in Fig. 5e exhibits a peak around 400 cm −1 , a shoulder around 100 cm −1 and a slow decay for lower frequencies. The peak around 400 cm −1 (12 THz) is caused by oscillations of the oxygen-oxygen separation, R OO , which couple to the proton position d via the most likely proton position d * [R OO (t)]; in simple terms, the proton vibrates with the water molecule it is bound to. The oxygen vibrational time scale τ R OO = 86 fs, indicated in Fig. 5b, follows from the peak of the power spectrum of R OO around 400 cm −1 , which is plotted in the lower panel of Fig. 5e and agrees perfectly with the peak in ωχ 00 TW . This peak is in fact also well visible in our experimental THz/FTIR difference spectra, shown again in Fig. 6a as a broken red line for a 6 M HCl solution, the dotted red line shows the corresponding simulated difference spectrum. Note that this translational vibration of two water oxygens in the transient H 5 O 2 + complex is about twice as fast as the translational vibration of two hydrogen-bonded water molecules in pure water, which gives rise to the well-known IR signature around 200 cm −1 , shown as a blue solid line in Fig. 6b obtained from pure-water simulations 12 . This frequency shift is the reason why the water-vibration peak appears prominently in the difference spectra in Fig. 6a.
The shoulder in ωχ 00 TW around 100 cm −1 is related to the transfer waiting time τ TW , which is the average time between two consecutive complete proton-transfer events, as predicted from an analytically solvable barrier-crossing model 38 . In Fig. 5f we show distributions of the transfer-waiting first-passage times, i.e., distributions of the time difference between crossing the most likely proton position d * [R OO (t)] on one side of the midplane d = 0 and crossing d * [R OO (t)] on the other side of the midplane d = 0 for the first time, for the three HCl concentrations. The distributions are essentially exponential in nature, which means that transfer events occur at a roughly constant rate and reflects the stochastic nature of the excess-proton transfer process in this reduced one-dimensional description. The mean of these firstpassage distributions defines the transfer-waiting time τ TW , which is given in Table 1 and increases with rising HCl concentration. This indicates that hydronium ions have a slightly longer life time at higher HCl concentrations. In contrast, both complete and incomplete transfer-path times interestingly show no dependence on the HCl concentration. The inverse of the transfer-waiting time τ TW = 280 fs for 6 M, which is shown in Fig. 5e as a vertical line, is located at 120 cm −1 (3.5 THz) and corresponds well to the position of the shoulder in ωχ 00 TW , which confirms the connection between the transfer-waiting time and the spectroscopic signature around 100 cm −1 that is predicted by analytical theory 38 . We note that the total length of the continuous proton-transfer trajectories are roughly twice as long as the mean transfer-waiting times, meaning that typically a few back-and-forth protontransfer events occur in each trajectory (see Supplementary Note 7 for more details).
The characteristic time scales of each contribution, i.e., the transfer-waiting time τ TW = 280 fs, the water-oxygen vibrational time τ R OO = 86 fs, the transfer-path times τ c TP = 14.1 fs and τ i TP = 12.7 fs and the normal-mode times τ NM = 11-17 fs are unambiguously extracted from the simulations and characterize both the trajectory contributions in the time domain in Fig. 5a-d, where they are included as horizontal black bars, and also the different spectral contributions in Fig. 5e.
We comment on the subtle spectral features in the range 1400-1800 cm −1 in Fig. 5e, where small but distinct peaks are revealed in the different spectral contributions. The transferwaiting contribution (blue line) peaks at about 1650 cm −1 , hinting to a weak coupling to an unperturbed water bending mode of the flanking water molecules. The transfer-path contribution (red line) peaks at 1750 cm −1 , the location of the experimental acid bend signature, which suggests that the acid bend couples particularly to the transfer path motion of the excess proton. Note, that even though the acid bend is primarily produced by the excess-proton motion orthogonal to the d coordinate, this contribution to the isotropic spectrum is largely compensated by motion of the flanking water molecules, as shown in Supplementary Fig. 11. The normal-mode contribution (green line) peaks around 1500 cm −1 , consistent with previously calculated normal-mode spectra of Eigen-like solvated proton structures 25 . We thus see that our trajectory-decomposition technique also allows to disentangle the various normal-modes obtained for the solvated excess-proton complex.
So far we have concentrated on the excess-proton spectral contribution and not discussed the chloride contribution. The decomposition of the total simulated 6 M HCl spectrum in Fig. 6b (red solid line) into the chloride contribution (green solid line, including all cross correlations) and the remainder (gray solid line) demonstrates a prominent chloride peak around 150 cm −1 , which translates to a corresponding time scale of τ Cl À = 185 ps and is due to the rattling of a chloride in its hydration cage 9,10 . This peak is also seen in the simulated difference spectrum in Fig. 6a (red dotted line) and is slightly shifted to larger frequencies in the experimental difference spectrum (red broken line). At 100-200 cm −1 the remainder contribution in Fig. 6b (gray solid line) is significantly stronger than the pure-water spectrum (blue line), indicating that a process related to excess-proton motion significantly contributes in this wavenumber range. We suggest that this process is the excessproton transfer-waiting contribution, which in Fig. 5e is shown to produce a broad shoulder around 100 cm −1 .
While the generally good agreement between our simulated and experimental spectra supports our chosen simulation methodology, it is clear that our classical treatment of nuclei motion is a drastic approximation and therefore some of the agreement might be due to fortuitous cancellation of errors. Interestingly, previous studies found no significant differences between IR spectra computed from simulations with and without nuclear quantum effects below 3000 cm − 1 34,35,55 , which might suggest that quantum-mechanical zero-point motion influences the excess-proton dynamics less than the instantaneous excessproton distribution. We discuss quantum nuclear effects and basis set issues in the "Methods" section and Supplementary Note 1. Furthermore, we also compare our simulation results for radial distribution functions 28,33,52,56,57 (Supplementary Note 4) and proton diffusion coefficients 28,33,[58][59][60]

Discussion
We show that the spectroscopic signature of proton-transfer dynamics between two water molecules in hydrochloric acid (HCl) solutions can be investigated by trajectory decomposition into transfer-waiting (characterized by the time scale τ TW ), transfer-path (characterized by τ TP ) and normal-mode contributions (characterized by τ NM ). The decomposition is performed in the two-dimensional coordinate system that is spanned by the excess-proton position and the oxygen-oxygen distance of the two flanking water molecules and operates both in the time domain as well as in the frequency domain. The coupling of the excess-proton motion to the relative oscillations of the two flanking water molecules produces a fourth spectral protondynamics contribution (characterized by τ R OO ). The dynamics of each of the four contributions are described by distinct time scales with the ordering τ TW > τ R OO > τ TP > τ NM and therefore contribute with distinct peaks to the excess-proton IR spectrum. Due to the high proton charge, the excess-proton spectrum, and due to correlation effects with the neighboring water molecules in particular the excess-proton motion along the axis connecting the oxygens, contributes significantly to the IR difference spectrum of HCl. Our experimental THz/FTIR difference spectra resolve the slowest time scales, τ TW and τ R OO , of which the former one is overlaid by an additional spectral contribution due to rattling of the chloride ions, characterized by yet another time scale τ Cl À , which is close to τ TW . Mid-IR experimental difference spectra from literature on the other hand are compatible with our predicted spectra associated with the τ TP and τ NM time scales.
In contrast to the transfer-path time τ TP , the transfer-waiting time τ TW shows a weak dependence on the HCl concentration. Possible reasons include ionic screening and repulsion effects between neighboring excess protons, but also entropic effects due to the reduced number of accepting water molecules have been discussed 28,32,62 . This is consistent with experimental results showing a decrease of the excess-proton diffusivity with increasing HCl concentration 60 , which is reproduced in our simulations (see Supplementary Fig. 19). One should note that the transfer-waiting times in local transient H 5 O 2 + complexes include back-and-forth proton-transfer events, which do not contribute to the long-time excess-proton diffusion 28,33,58,59,63 but nevertheless have a pronounced spectroscopic signature 38 .
Our results nicely complement recent normal-mode calculations. We find that the continuum band stems from normal-  The standard errors of the transfer-waiting time τ TW and the transfer-path times τ i TP and τ c TP are estimated from the variances of the fitted distributions. The errors of the normal-mode times τ NM , the oscillation time of the two flanking water molecules τ R OO and the rattling time of chloride ions τ Cl À are estimated from the resolution of the underlying spectra. mode vibrations of less symmetric, i.e., more Eigen-like, configurations of the excess proton, which are strongly influenced by the oxygen-oxygen separation R OO 14,19,25,26 . On the other hand, our transfer-path signature, which is dominated by a broad absorption around 1200 cm −1 , shows striking similarity with normal-mode spectra computed for more symmetric, i.e., more Zundel-like, configurations of the excess proton 14,19,25,26 . It is not implausible that normal modes for small separations of the two flanking water molecules show similar spectroscopic signatures as the transfer paths we extract from our simulation trajectories. Yet, in a normal-mode picture the interconversion between the metastable Zundel-like and Eigen-like excess-proton states cannot be explained consistently, even though the importance of this process for a complete description of the mid-IR signatures was acknowledged several times 14,[24][25][26]35 . In fact, the broad distribution of this interconversion time scale is demonstrated by the transfer-waiting time distribution and also gives rise to a distinct spectral signature, that we identify in the THz regime.
A recent study employing similar simulation techniques has decomposed the proton power spectra with respect to the proton asymmetry coordinate and thereby reached similar conclusions to ours: Eigen-like configurations give rise to the continuum band while Zundel-like configurations dominantly contribute around 1200 cm −1 47 . That study also determined the proton-transfer time scale using two-dimensional transition state theory and Marcus theory of ion pairing and finds this time scale to be concentration dependent, in agreement with our and previous observations.
In summary, in many theoretical treatments, only normal modes of meta-stable or stable states are assumed to produce spectral contributions. Any spectral mode is therefore interpreted as being due to a meta-stable state, consequently, broad spectral modes are often interpreted as reflecting a wide collection of normal modes with slightly different frequencies. In this paper, we show that transfer and barrier-crossing events of charged particles as well as transfer paths create spectral features by a mechanism that is very different from a normal-mode picture and that these spectral features are broadened by the stochastic nature of the transfer dynamics 38 . In fact, the strength or frequency of a spectral feature does not allow to tell whether it is caused by normal-mode oscillations in a stable or meta-stable state or whether it is caused by transfer or barrier-crossing dynamics.

Methods
Computational methods: ab initio molecular dynamics simulations. The Born-Oppenheimer ab initio MD simulations of pure water and HCl solutions at three different concentrations were performed with the CP2K 7.1 software package using a polarizable double-zeta basis set for the valence electrons, optimized for small molecules and short ranges (DZVP-MOLOPT-SR-GTH, with the exception of the chloride anions, that were modeled including diffuse functions in the aug-DZVP-GTH basis set), dual-space pseudopotentials, the BLYP exchangecorrelation functional, and D3 dispersion corrections [64][65][66][67][68] . The cutoff for the plane-wave representation was 400 Ry. The system parameters are summarized in Table 2.
Before production, each system was equilibrated in classical MD simulations for 200 ps under NPT conditions at atmospheric pressure and 800 ps under NVT conditions at 300 K, using the GROMACS 2020.5 software 69 with the SPC/E water model 70 . The force fields for Cl − and H 3 O + were taken from 71 . The ab initio MD simulations were subsequently performed using a time step of 0.5 fs under NVT conditions at 300 K by coupling the system to a CSVR thermostat with a time constant of 100 fs 72 .
Dipole moments were obtained after Wannier-center localization of the electron density at a time resolution of 2 fs or 4 fs. At each time step, the Wannier centers were assigned to the closest oxygen or chloride ion. Water molecules were assembled by assigning each proton to the closest oxygen nucleus, thereby forming either water or hydronium ions. For the hydronium ions, all protons were treated as excess-proton candidates and further processed based on a dynamical criterion as discussed in the main text and Supplementary Methods 2. The dipole moments p follow as a sum over the respective position vectors r i and charges q i (q = 2 e for Wannier centers, and reduced core charges for nuclei), p = ∑ i q i r i for the whole or desired sub systems.
Linear response theory relates the dielectric susceptibility χ(t) to the equilibrium autocorrelation of the dipole moment C(t) = ∑ D 〈p(t)p(0)〉, reading in Fourier space with system volume V, thermal energy k B T, vacuum permittivity ϵ 0 and D being the number of Cartesian dimensions of the polarization vector p. IR spectra can therefore be calculated straight-forwardly from sufficiently sampled trajectories of the ab initio MD simulation data using Eq. (4)  wherepðωÞ is the Fourier-transformed dipole-moment trajectory with length L t and the asterisk denotes the complex conjugate. Alternatively, for charged subsystems, as in the case of the chloride ions, the computation using the time derivative of the polarization, i.e., the current j ¼ d dt pðtÞ, is preferable Quantum corrections have previously been addressed 73 , but were not applied here.
Since the Wannier-center localization time step Δt WC = 4 fs is larger than the original simulation time step Δt = 0.5 fs, the analysis is performed on two types of trajectories stemming from the same simulations: one set of trajectories containing the electronic degrees of freedom and another set of trajectories of higher time resolution but only containing nuclei positions. This higher resolution data is used for the calculation of excess-proton spectra and kinetics.
All spectra were smoothed by a convolution with a Gaussian kernel of varying width, depending on their respective resolution. We used a standard deviation of 55, 20, and 50 cm −1 for bulk spectra, the H 5 O 2 + complex difference spectrum in Fig. 3c and excess-proton spectra, respectively. Experimental data was smoothed using a standard deviation of 3 cm −1 .
To address the quality of the chosen basis set, shorter simulations at 6 M were performed using the non-short range basis set (DZVP-MOLOPT-GTH) as well as a triple-zeta doubly polarizable (TZV2P-GTH) basis set. Spatial correlations in the data are compared in Supplementary Fig. 9. While the coordination of excess protons with chloride ions slightly increases when the more elaborate basis sets are used, no significant differences in correlations between excess protons and oxygen nuclei are found, which are the focus of this study.
Experimental methods: THz absorption measurements. THz spectroscopic measurements in the 30-650 cm −1 frequency range were done with a commercial Fourier Transform spectrometer (Bruker Vertex 80v, Germany) equipped with a mercury light source and a liquid helium cooled bolometer detector (Infrared Laboratories, Germany). Spectra result from an average of 128 scans with a resolution of 2 cm −1 . The liquid sample cell is composed of diamond windows (Diamond Materials GmbH, Germany) in which a Kapton spacer of approximately 13 μm was placed between the windows to fix the sample thickness. The exact thickness of the sample cell was determined from the etaloning pattern of the empty sample cell. The temperature of the sample was held constant at 20.0 ± 0. 2 ∘ C by an external chiller. The measured frequency-dependent extinction coefficient, α solution (ω), is determined using the Beer-Lamber law where d is the sample thickness, I water (ω) and I solution (ω) are the experimental transmitted intensities of the water reference and the sample. α water (ω) is the extinction coefficient of bulk water and is taken from literature 74 . The extinction coefficient α(ω) is converted to the absorption spectrum, proportional to the imaginary part of the dielectric susceptibility χ″(ω), by fitting the spectra and performing a Kramers-Kronig transform as presented in Supplementary Methods 1. The difference absorption spectra of HCl solutions with respect to pure water and normalized with respect to the water concentration are given by where c w and c 0 w are the concentration of water in the aqueous HCl solutions and bulk water, respectively, determined from the solution density at room temperature.

Data availability
The datasets generated and analyzed during the current study are available from the corresponding author on request.